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ABSTRACT 

T-REx (Tau Retrieval of Exoplanets) is a novel, fully Bayesian atmospheric retrieval code cus¬ 
tom built for extrasolar atmospheres. In Waldmann et al. (2015) the transmission spectroscopic 
case was introduced, here we present the emission spectroscopy spectral retrieval for the T-REx 
framework. Compared to transmission spectroscopy, the emission case is often significantly more 
degenerate due to the need to retrieve the full atmospheric temperature-pressure (TP) profile. 
This is particularly true in the case of current measurements of exoplanetary atmospheres, which 
are either of low signal-to-noise, low spectral resolution or both. Here we present a new way of 
combining two existing approaches to the modelling of the said TP profile: 1) the parametric 
profile, where the atmospheric TP structure is analytically approximated by a few model param¬ 
eters, 2) the Layer-by-Layer approach, where individual atmospheric layers are modelled. Both 
these approaches have distinct advantages and disadvantages in terms of convergence properties 
and potential model biases. The T-REx hybrid model presented here is a new two-stage TP 
profile retrieval, which combines the robustness of the analytic solution with the accuracy of 
the Layer-by-Layer approach. The retrieval process is demonstrated using simulations of the 
hot-Jupiter WASP-76b and the hot SuperEarth 55 Cnc e, as well as on the secondary eclipse 
measurements of HD 189733b. 


Subject headings: methods: data analysis — methods: statistical — techniques: spectroscopic — radia¬ 
tive transfer 


1. Introduction 

The characterisation of extrasolar planets 
through the spectroscopic measurements of their 
atmospheres has become a well established field 
today (Tinetti et al. 2013). In Waldmann et al. 
(2015), we presented the T-REx (Tau Retrieval of 
Exoplanets) suite for transmission spectroscopic 
measurements. In this paper we introduce the 
corresponding retrieval for emission spectroscopic 
data. 

Transmission and emission spectroscopy carry 
highly complementary aspects. Whereas trans¬ 
mission spectroscopy is less sensitive to thermal 
gradients, the emission spectroscopy case probes 
the full temperature-pressure profile (TP-profile 
hereafter) of the atmosphere. This makes the 
emission case significantly harder to constrain 
without the luxury of in-situ measurements. King 


(1958) was the first to suggest remote sensing 
of the planetary atmospheric temperature struc¬ 
ture through the infra-red, thermal radiance of 
the planet. Kaplan (1959) expanded on this pio¬ 
neering work by laying the foundations of retriev¬ 
ing exact TP-profiles from emission spectroscopic 
measurements. Remote sensing of planetary at¬ 
mospheres in our solar system has been a long 
story of success (e.g. Wark & Hilleary 1969; Con- 
rath et al. 1970; Hanel et al. 1972; Conrath et al. 
1973; Rodgers 1976; Hanel et al. 1981; Eletcher 
et al. 2007; Irwin et al. 2008) as well as Hanel et al. 
(2003), and references therein. More recently, 
emission spectroscopy remote sensing has been ex¬ 
panded to exoplanetary atmospheres (Madhusud- 
han & Seager 2009; Lee et al. 2011; Line et al. 
2012; Griffith 2014). Line et al. (2013) provides 
a comparison of existing exoplanetary emission 
spectroscopy retrieval codes. 
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1.1. r-REx 

T-REx is a novel, Bayesian atmospheric spec¬ 
troscopy retrieval suite designed for extrasolar 
planets. In Waldmann et al. (2015), hereafter 
W15, we introduce the overall architecture, data 
handling, minimisation/sampling routines, han¬ 
dling of molecular line lists, Bayesian model se¬ 
lection and the transmission spectroscopy forward 
model. In this publication we present the emis¬ 
sion spectroscopy forward model as well as the 
atmospheric temperature-pressure (TP) profile re¬ 
trieval implemented. 

Figure 1 shows a schematic diagram of the T- 
REx architecture for the emission retrieval. The 
program can be subdivided into five sections: 

1. Inputs: Program inputs such as parame¬ 
ter files, molecular line-lists (ExoMol Ten¬ 
nyson & Yurchenko 2012), (HITEMP Roth¬ 
man et al. 2010), stellar spectra (PHOENIX, 
Allard et al. 2012) and spectroscopic obser¬ 
vations to be analysed. 

2. Model & Data handling: The Central Data 
Module manages all calls to the TP-Profile 
and Forward Model Modules and provides a 
standardised, common interface to different 
sampling routines implemented. 

3. Sampling: T-REx features three indepen¬ 
dent sampling routines: Maximum Likeli¬ 
hood estimation (MLE, based on the LM- 
BFGS minimisation in W15), Markov Chain 
Monte Carlo (MCMC) and Nested Sampling 
(NS) routines. We refer the reader to W15 
for implementation details. 

4. Post analysis & refinement: This stage dif¬ 
fers from W15 by providing a two itera¬ 
tion stages for the determination of the TP- 
Profile, see section 3. 



Output 


Fig. 1.— Flowchart illustrating the modular de¬ 
sign of T-REx. This flowchart is an update from 
Waldmann et al. (2015), reflecting the differing ar¬ 
chitecture of the emission spectroscopy code. 

Chandrasekar 1960; Coody & Yung 1989; Lion 
2002; Hanel et al. 2003; Mihalas & Mihalas 2013). 
In what follows, we closely follow the nomencla¬ 
ture of Lion (2002). For non-scattering atmo¬ 
spheres in local thermodynamic equilibrium, the 
basic equation describing the thermal radiation 
of an atmosphere is given by the Schwartzschild 
equation: 

( 1 ) 

where T is the intensity per wavelength. A, Bx is 
the Planck function at temperature T, /i = cos^ 
is the upwards inclination, r is the overall optical 
depth, given as function of altitude (z) by 


5. Output: The final posterior distributions are 
analysed and the final model spectrum, TP- 
Profile and mixing ratios are returned. 

2. Forward Model 

We briefly describe the radiative transfer for¬ 
ward model, for a more exhaustive discussion we 
refer the reader to the standard literature (e.g.. 


= Y (2) 

m=l 

where rA,m denotes the optical depth per absorb¬ 
ing species, m, given by 

r^oo 

T\m = / •iX,m{z')Xm{z')pN{z')<^z' . (3) 

J z 

Here is the absorption cross section, Xm Ihe 
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column density and pN the number density. We 
can now express the upwards welling radiance as 


Jr 


( 4 ) 

where the first right-hand-side term is the radi¬ 
ation at the planetary surface (or defined surface 
pressure for gaseous planets), and the second term 
denotes the integrated emission contributions for 
individual plane-parallel layers. The monochro¬ 
matic transmittance and its derivative (weighting 
function) can be defined as 


T\{t) = e 


dTxir) 

dr 


—e 


— T 


( 5 ) 


Hence we express the total integrated radiation at 
the top of the atmosphere (TOA, r = 0, z = oo) 
as 


Ix{t = 0 ) = Bx{T,)e--^ + Bx{Tr)^^^^dT 

(6) 

where and Tg are the optical depth and temper¬ 
ature at the planetary surface. The final exoplan¬ 
etary emission spectrum is now given by 


Fp/F, = 


Ix{t = 0 ) 
h 


X 



( 7 ) 


here F is the stellar intensity. By default T - 
REx interpolates the stellar flux from a library of 
Phoenix^ stellar spectra, gridded at lOOK inter¬ 
vals. Alternatively T-REx can approximate the 
stellar intensity using a black-body at the stellar 
temperature. 


3. Temperature - Pressure Profile 

The determination of the vertical atmospheric 
temperature profile (also referred to as temperature- 
pressure profile or TP-profile for short) is one of 
the key challenges in the retrieval of atmospheric 
emission spectra. Typically two approaches exist 
in the retrieval of the TP-profile: 1) Layer-by-layer 
retrieval, 2) Analytic parameterisation. 


^http://WWW.hs.uni-hamburg.de/EN/For/ThA/phoenix/ 
index.html 


1) The layer-by-layer approach consists of mod¬ 
elling the temperature of each plane-parallel at¬ 
mospheric layer independently. This approach is 
usually adopted for retrieval work of the Earth’s 
atmosphere and solar system planets (Hanel et al. 
2003; Rodgers 2000). The advantage of such an 
approach is its non-parametric nature, i.e. no 
prior knowledge is imposed on the temperature 
retrieval nor is the solution limited by potential 
restrictions and biases of an analytic TP-profile. 
The obvious disadvantage lies in the potentially 
poor convergence properties of such an approach 
in low S/N scenarios. Often, data quality or sparse 
spectral sampling do not allow for a pure layer- 
by-layer retrieval due to the large number of free 
parameters incurred. This is particularly true for 
current exoplanet spectroscopy data, which is ei¬ 
ther of relatively low S/N, low spectral resolution, 
or both. A common approach adopted is to impose 
a ‘regularisation’ of the TP-profile based on the 
physical reality that adjacent atmosphere layers 
should exhibit some correlation in temperature. 

2) The second approach to the TP-profile re¬ 
trieval is to adopt an analytic (here also referred 
to as ‘parametric’) model. Such models analyti¬ 
cally approximate the mean underlying physics of 
the atmospheric thermal structure. Such models 
have the advantage of containing far fewer free pa¬ 
rameters compared to the layer-by-layer approach, 
hence they converge faster. However, the solution 
is constrained within the bounds of the model as¬ 
sumed. 

In summary, the layer-by-layer methodology 
is most objective but features poor convergence 
properties, whilst the parametric model is less ob¬ 
jective but converges more robustly. Eor a review 
of existing implementations of both approaches in 
the field of exoplanet atmospheric retrieval we re¬ 
fer the reader to Line et al. (2013). 

As described in section 3.2. T-REx employs 
a hybrid model combining both, parametric and 
layer-by-layer approaches in a two stage retrieval 
process. In section 3.1 we briefly describe the 
parametric models implemented in T-REx as part 
of the T-REx retrieval framework. 

3.1. Parametric Models in the literature 

Analytical global-average TP-profiles exist in 
various flavours ranging from 2-stream purely 
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radiative to radiative-convective approximations 
and global circulation models (GCMs). We re¬ 
fer the reader to the relevant literature for a 
more in-depth discussion on the various ana¬ 
lytic approaches (e.g. Lion 2002 ; Hubeny et ah 
2003; Hansen 2008; Burrows et al. 2008; Show¬ 
man et al. 2009; Madhusudhan & Seager 2009; 
Pierrehumbert 2010; Guillot 2010; Heng et al. 
2012 , 2014; Robinson & Gatling 2012 ). For the 
Stage 1 approach, T-REx features a purely radia¬ 
tive 2 -stream model as well as a choice of simpler 
TP-profiles based on purely geometric considera¬ 
tions (see 3.1.1). As described in section 3.2, the 
exact form of the TP-profile is less relevant in our 
case as the results will be refined in a second stage 
fitting. 

Following Guillot (2010), the mean global tem¬ 
perature profile for a simple radiative downstream- 
upstream approximation can be expressed as 


T\t) 




+ 



(8) 


where Tint is the planet internal heat flux, Ti^r the 
solar fiux at the top of the atmosphere and 
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I.(f-l) 


-nr 


( 9 ) 




where 71 = is the ratio of mean opacities 

in the optical and infra-red {hZiR) and E 2 is 
the second-order exponential integral given by 

En+liz) = l[e“^ - zEn{z)]. (10) 

n 

We note that similar parameterisations exist in the 
literature (e.g., eq. (18), Robinson & Gatling 
2012 ; Pierrehumbert 2010). We also include the 
variation by Line et al. ( 2012 ) and Parmentier 
& Guillot (2014) including two optical opacity 
sources Rm and and a weighting factor be¬ 
tween optical opacities (left as free parameter) a. 


STL {2 


STL 


T\r)=T^[-+T\ + T^{l- a)L{r) 


4 V3 


OrpL 


( 11 ) 


The temperature as function of opacity r can be 
mapped to a pressure grid by assuming the follow¬ 
ing relation 


i^irP 
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( 12 ) 


3.1.1. Other TP-profiles 

In addition to the above TP-profiles, we include 
an isothermal profile as well as a ‘3-point’ and ‘4- 
point’ profile. The 3-point profile is purely geo¬ 
metric and keeps the top of atmosphere tempera¬ 
ture, Ttoa^ the tropopause temperature and pres¬ 
sure, Ti, Pi, and the surface (or 10 bar pressure) 
temperature Tio^ar as free variables. The TP- 
profile is then linearly interpolated in ln(P). The 
4-point profile adds an extra variable temperature- 
pressure point to the profile. 


3.2. The T-REx Hybrid Model 



Temperature 


Fig. 2 .— Left: A given Temperature - Pressure 
Profile without inversion. Right: Correlation ma¬ 
trix C, equation 14, of the TP-profile on the left. 
Atmospheric layers with similar temperatures fea¬ 
ture a high layer-layer correlation whilst layers 
with temperature differences feature a lower cor¬ 
relation. 

In T-REx we have worked towards a hybrid 
solution of the above mentioned approaches in 


4 


















which the final solution is not bound by a para¬ 
metric model but does not inherit the potentially 
poor convergence properties (and susceptibility to 
noise) of a layer-by-layer approach. Here we pro¬ 
ceed in two stages: 1) We perform a parametric 
model retrieval, 2) we take the retrieval solution 
to ‘guide’ a second layer-by-layer retrieval, relax¬ 
ing the parametric model constraint and thereby 
fine-tuning the TP-profile. 

3.2.1. Stage 1 

In Stage 1 we adopt a classical parametric 
model retrieval using one of the TP-profile param- 
eterisations implemented in T-REx. The idea is 
to constrain an initial best-fit of the model and 
whilst a realistic model (i.e. one well suited to de¬ 
scribe the underlying TP-profile) is preferable, it 
is not a hard requirement. We now fit the solution 
using either of T-REx’s sampling routines (MLE, 
MCMC, NS). The error on the sampled paramet¬ 
ric model parameters is converted to one a lower 
and upper temperature bounds for a layer-by-layer 
model. This is done using a numerical model per¬ 
mutation analysis of the Stage 1 parameters. 

We then calculate the following distance matrix 


3.2.2. Stage 2 

In the second stage of our TP-profile retrieval 
we relax the original solution found in Stage 1 and 
by that ‘fine-tune’ the TP-profile. We construct a 
second correlation matrix imposing an exponential 
correlation length across pressure levels (Rodgers 
1976, 2000) 

Sij = {SiiSjj)^^^eyip (15) 

where c is the correlation length in terms of at¬ 
mospheric scale heights. We set c = 3.0 by de¬ 
fault, but can otherwise be user defined. Rodgers 
(2000) gives an advisable range between 1.0 - 8.0, 
with Irwin et al. (2008) using a default of c = 1.5 
and Line et al. (2013) c = 7.0. Larger values of 
c correspond to a stronger regularisation of the 
TP-profile (i.e. smoothing). A stronger regular¬ 
isation can be useful in poorly constrained data 
sets (either low S/N, sparsely sampled or both). 
Equation 15 is identical to that used in the layer- 
by-layer approach by Irwin et al. (2008). We now 
construct a hybrid correlation matrix by combin¬ 
ing equations 14 & 15 


Al = \f,-f,\^ + {a, + aj)^ (13) 


'Dij{a) — aCij + (1 — OL)Sij (16) 


where T is the maximum likelihood temperature 
estimator of the parametric model fit for the 
and atmospheric layer and a is the respective 
error. Equation 13 captures the layer to layer tem¬ 
perature variations in the TP-profile and is hence 
conceptually similar to the Jacobian of the pro¬ 
file. We normalise in terms of the minimal 
and maximal temperature variations found in the 
TP-profile, 


^ Ain — argmin A 

Cij = 1.0 - — 

argmax A 


(14) 


which can be understood as a positively defined 
temperature correlation matrix with layers most 
similar in temperature featuring the highest cor¬ 
relation and layers most dissimilar resulting in a 
very low correlation. An example of C for a given 
TP-profile can be found in fig. 2. 


where a is a scaling factor ranging from zero to 
one. Eigure 3 shows instances of V{a) at varying 
values of a. We now correlate our layer-by-layer 
TP-profile with 'D(a) whereby leaving a as free pa¬ 
rameter to be fitted. By allowing a to vary, we can 
dynamically relax the parametric model solution 
from Stage i, a = 1.0, to an unconstrained solu¬ 
tion, Of = 0.0. The advantage of such an approach 
is twofold: 1) Since the Stage 1 fitting should have 
already achieved a solution close to the real max¬ 
imum likelihood, convergence in the second stage 
towards the volume of lowest is significantly 
faster; 2) T-REx can freely decide to relax the 
Stage 1 solution should this be favoured by the 
data. Practically this happens frequently at later 
stages of the fitting. Whereas T-REx initially con¬ 
verges towards the Stage 1 solution (i.e. a ^ 1), 
at later stages of the fitting the code begins to 
reject the parametric model (i.e. a ^ 0) as it 
‘fine-tunes’ the original solution. This ‘tuning’ can 
achieve significantly higher maximum likelihoods 
than the Stage 1 fitting alone. 
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Fig. 3.— Hybrid correlation matrix P, equation 16, at different values of a. The left most is the pure 
Stage 1 correlation matrix, C, whereas the right plot is the pure correlation-length-only matrix S', eq. (15). 
Over-plotted are contours of V{a). The middle panel shows an intermediate state of left (60%) and right 
(40%). 


3.3. Non-linear Sampling 

In the approach explained in section 3.2, we 
keep an even sampling of atmospheric layers in 
log(P) for Stage 2. For well sampled and high 
S/N data, this approach is adequate. However for 
coarsely sampled and/or poor S/N data, it is of¬ 
ten advisable to reduce the number of free param¬ 
eters to a minimum to aid convergence. In these 
cases we can utilise the sparsity of the Stage 1 
TP-profile solution to devise a nonlinearly sam¬ 
pling of the exoplanetary atmosphere. We base 
our compression algorithm on the fact that only 
changes in the temperature-pressure gradient need 
to be modelled, i.e., an isothermal TP-profile can 
be perfectly described by a single temperature pa¬ 
rameter whereas areas of non-isothermal structure 
need more parameters to capture variation. The 
compression algorithm uses the Stage 1 correla¬ 
tion matrix C to only retain TP-profile layers cor¬ 
responding to a > 2% change in the TP gradi¬ 
ent with respect to the previously retained layer. 
We graphically depict this in figure 4. In addition 
to TP-profile layers sampled this way, we also in¬ 
clude an extra sampling layer whenever no change 
in thermal gradient was detected for > 10 layers. 
The majority of the TP-profile variations should 
be captured by the Stage I retrieval and Stage 2 
is ‘fine-tuning’ this solution. The inclusion of a 
coarse sampling in isothermal regimes does allows 
the Stage 2 retrieval to deviate from any Stage I 



Fig. 4.— Nonlinear TP-profile sampling on cor¬ 
relation matrix C (same as right hand plot in fig¬ 
ure 2). Starting at the top of the atmosphere (red 
dot) we retain all layers in the TP-profile that cor¬ 
respond to a change in gradient > 2% with respect 
to the previously retained layer until the bottom 
layer (green dot) is reached. The selected layers 
are denoted by white dots and arrows represent 
the path of the compression algorithm across the 
correlation matrix C. Should no gradient change 
be detected for > 10 layers, an extra sampling 
point is introduced (black dots). 
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solution, should this be supported by a higher 
global likelihood. Using the non-linear sampling, 
we can reduce a 100 layer atmospheric model to 
typically 15-25 free parameters. 

4. Validation of T-REx 

We demonstrate the emission spectroscopy re¬ 
trieval with two simulations and the secondary 
eclipse spectrum of HD189733b. The simulations 
are as follows: 1) high S/N observation simulation 
of a hot Jupiter similar to WASP-76b; 2) observa¬ 
tions of the hot SuperEarth 55 Cnc e, simulated 
for a 1 meter class spectroscopic space mission. 
In our simulations we opt for two oxidised atmo¬ 
spheres at high temperatures (> 1500K). 

For each retrieval stage we calculate the global 
Bayesian Evidence of the solution set. Here the 
Bayesian Evidence (E, or simply Evidence here¬ 
after) is given as the integral of the product of the 
global likelihood and the prior space 

E = J p{e\M)P{x\e,M)de (i?) 

where P{6\Ai) is the prior over the parameter 
space 0, for retrieval model M. P(x|0) is the like¬ 
lihood for the observed data vector x given the 
parameter space and retrieval model. Retrieval 
Evidences are reported in tables 1 and 2. 

4.1. WASP-76b 

WASP-76b (West et al. 2013) is a very hot- 
Jupiter orbiting a late F7 star. It is highly in¬ 
flated at 1.83lo‘.o4 ^Jupi 0-92 ± 0.03 Mjup and 
Tequ ^ 2200K. We take its bulk and orbital prop¬ 
erties and generate a simulated observation at a 
resolution of 300, spectral range of 0.5 - 20/im and 
lOOppm errors. We set the atmospheric composi¬ 
tion to 1 X 10“"^ H 2 O, 1 X 10“^ CO and 1 x 10“^ 
CO 2 . The input TP-profile is shown in figure 7 
(black line). It is important to note that here (as 
well as in section 4.2) the input TP-profile was 
generated using a script external to T-REx with 
no relation to either Stage 1 parametrisation. This 
provides an adequate test-bed for the Stage 1 
fitting to accurately retrieve an arbitrary atmo¬ 
spheric profile. 

As described in section 3.2, we retrieve the TP- 
profile and abundances in two stages. For compu¬ 
tational efficiency reasons, here we only compute 



0 10 20 30 40 50 60 70 


Layers 

Fig. 5.— Correlation matrix, C, derived from 
Stage 1 TP-profile shown in figure 7. The in¬ 
put model TP-profile is shown as black continuous 
line. Otherwise identical to figure 3. 



Fig. 6.— showing the input spectrum of a WASP- 
76b type hot-Jupiter, gray, and Stage 1 & 2 fit¬ 
ting on the top and bottom respectively. Both 
fits converged but the Stage 2 fit reached a higher 
maximum likelihood. 


7 










Table 1: Retrieved abundances for hot Jupiter WASP-76b. Top row: log Bayesian Evidence for 
Stage 1 & Stage 2 models. Differences above |Alog(E) = 5| are very significant. Here Alog(E) = +80.12 
indicating a significantly improved fit in Stage 2. 



Input Model 

Stage I Retrieval 

Stage 2 Retrieval 

log{E) 

NA 

-43.92 

36.20 

log(H20) 

-4.0 

-3.458 ± 0.104 

-3.660 ± 0.107 

log(CO) 

-4.69 

-4.352 ± 0.143 

-4.548 ± 0.II3 

log(C02) 

-6.0 

-5.374 ± 0.120 

-5.428 ± 0.II4 



Fig. 7.— TP-profiles for Stage 1 k. 2 results 
for WASP-76b in figure 6. Solid lines represent 
the mean and shaded regions the one sigma error 
bars. Solid black lines are the input TP-profile. 
Stage 2 takes the initial parametric TP-profile fit 
of Stage 1 and relaxes the solution. Four solu¬ 
tions were obtained for Stage 1, the highest max¬ 
imum likelihood solution (blue) traces the input 
TP-profile well, whilst local maxima underesti¬ 
mate the bulk temperature, see text. The Stage 2 
solution feature a significantly lower (or higher 
global Evidence) than all Stage 1 solutions. 


the Nested Sampling solutions (which are also the 
most accurate, see W15). Tests were run with 
both maximum-likelihood retrievals and MCMC 
retrievals and solutions between all sets of solu¬ 
tions are in good agreement. 

The Nested Sampling Stage 1 solution returns 
four local likelihood maxima of which the global 
maximum was selected. Figures 6 (top), 7 (left) 
and table 1 show the retrieved spectrum, TP- 
profile and retrieved abundances respectively. The 
Stage 1 TP-profile mostly captures the input TP- 
profile but shows unrealistic bumps and wiggles as 
well as unrealistic distributions of the one sigma 
error bar. These are artefacts of the parameter- 
isation in section 3.1. The three local maxima 
TP-profiles shown in figure 7 underestimate the 
bulk temperature of the planet, driving the re¬ 
trieval to assume unrealistically low abundances 
of molecular trace gases. In this example this is 
found to be a numerical effect that disappears by 
increasing the sampling grid density of the Nested 
Sampling. Nonetheless the potential degeneracy 
between TP-profile and molecular abundances is 
in atmospheric retrievals is worth noting. 

As described earlier, we computed the TP- 
profile covariance, C (figure 5), and tuned the 
Stage 1 retrieval by relaxing the parametric so¬ 
lution. Figures 6 (bottom), 7 (right) and table 1 
show the Stage 2 retrieval results. Inspecting fig¬ 
ure 6, both Stage 1 and Stage 2 retrievals fit the 
input spectrum sufficiently well but Stage 2 comes 
significantly closer to the ‘true’ TP-profile and 
trace gas abundances. Figure 8 shows the nor¬ 
malised emission contribution functions of both 
retrievals and their difference. This shows the 
planetary emission mainly emanating in the non- 
isothermal regions (up to the tropopause) of the 
TP-profile, as expected. However Stage 2 emis¬ 
sions originate significantly lower (by nearly an 
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order of magnitude) in the atmosphere (blue band 
in the bottom panel). 

4.2. 55 Cnc e 

We simulated an emission spectrum of the hot 
SuperEarth 55 Cnc e (Fischer et al. 2008). We use 
trace gas compositions of 1 x 10“"^ H 2 O, lx 10“^ 
CO and 1 x 10“^ CO 2 and a sharp TP-profile 
shown in fig. 10. As the previous WASP-76 ex¬ 
ample is currently unrealistically optimistic, given 
the combination of high S/N, moderate resolution 
(R r\j 300) and a very broad wavelength coverage, 
we opt for a more realistic example here. We calcu¬ 
lated the expected resolution and S/N for 100 co¬ 
added eclipses obtained by a one-meter-class tran¬ 
siting spectroscopy space-mission (e.g. the ESA 
M4 proposal ARIEL) over a wavelength range 
spanning 2 to 8/im. Using EChOSim^ an end- 
to-end simulator for transit spectroscopy space- 
missions (Pascale et al. 2014; Waldmann & Pascale 
2014) developed for the ESA M3 EChO proposal 
(Tinetti et al. 2012), we calculated realistic error 
bars for this hot SuperEarth, shown in figure 9. 

Similarly to section 4.1, figs 9, 10 & 11 show 
the Stage 1 & 2 retrieved spectra, TP-profiles and 
Stage 1 TP-profile covariance respectively. Table 2 
shows that Stage 2 retrieval converges at a sig¬ 
nificantly higher global Evidence and presents an 
improvement in the accuracy of abundances re¬ 
trieved as well as TP-profile retrieved. Figure 12 
shows the absolute difference between the Stage 1 
and Stage 2 model fits at a spectral resolution of 
1000. Here the discrepancies between both spec¬ 
tral fits are of the order of 2 x 10“^ or less. This 
is not very significant in terms of a naive fit 
to relatively low S/N data, but does significantly 
drive the retrieval of the TP-profile and trace gas 
abundances. This is reflected in the Bayesian Ev¬ 
idence measured for Stage 1 and Stage 2 mod¬ 
els. Figure 13 shows the contribution functions for 
both retrievals. Here we have positive and nega¬ 
tive temperature differences (see bottom left plot), 
resulting in a more complex contribution differ¬ 
ential (bottom left plot) than for WASP-76b. It 
should be noted that this contribution differen¬ 
tial is highly wavelength dependent and illustrates 
the sensitivity of varying wavelength ranges on the 
TP-profile retrievability. 


A 

stage 1 

HI.. ' ^ 

1^ 

stage 2 

. . 

... ' ^ 


stage 1 - Stage 2 I 

I 

.00 1.49 2.22 3.31 4.93 7.35 10.96 16.34 


Temperature (K) Wavelength (/im) 


Fig. 8.— Left, retrieved TP-profile, right: contri¬ 
bution functions for WASP-76b. First and second 
row: Best-fit TP-profiles and corresponding emis¬ 
sion contribution functions as function of pres¬ 
sure and wavelength for Stage 1 and Stage 2 re¬ 
spectively. Bottom row: Difference between both 
stages. Contribution difference is given as nor¬ 
malised fraction. 



Fig. 9.— Input spectrum of a 55 Cnc e type at¬ 
mosphere, gray, and Stage 1 & 2 fitting on the top 
and bottom respectively. Both fits converged but 
the Stage 2 fit reached a higher global Evidence. 
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Table 2: Retrieved abundances for hot SuperEarth 55 Cue e. Top row: log Bayesian Evidence for 
Stage 1 & Stage 2 models. 



Input Model 

Stage 1 Retrieval 

Stage 2 Retrieval 

log{E) 

NA 

75.40 

168.90 

log(H20) 

-4.0 

-4.168 ± 0.795 

-4.055 ± 0.571 

log(CO) 

-5.0 

-5.764 ± 1.248 

-5.613 ± 1.172 

log(C02) 

-5.0 

-5.236 ± 1.112 

-5.136 ± 1.019 



Eig. 10.— TP-profiles for Stage i & 2 results 
for 55 Cue e in figure 9. Both solutions converge 
within the calculated error bar. Stage 2 features 
a significant improvement in maximum likelihood 
achieved. 


0.205r 



Wavelength (/xm) 


Eig. 12.— showing the model fit difference be¬ 
tween Stage 1 and Stage 2 at spectral resolution 
of 1000. Model fit difference are of the order of 
2 X 10“^ or less. 




Eig. 13.— Contribution functions for 55Cnc e, 
otherwise same as figure 8. 


Eig. 11.— Correlation matrix, C, derived from 
Stage 1 TP-profile shown in figure 10. Otherwise 
identical to figure 5. 
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4.3. HD189733b 


Finally we test T-REx on the emission spec¬ 
trum of the hot Jupiter HD 189733b. Its emis¬ 
sion spectrum is amongst the most complete and 
best studied (e.g. Swain et ah 2009; Madhusud- 
han & Seager 2009; Lee et ah 2011; Line et ah 
2012 ; Line et ah 2014). We have compiled data 
sets ranging from the NIR (1.4/im) to the MIR 
(24/im), namely: HST/NICMOS (Swain et al. 
2009), Spitzer/IRAC 3.6, 4.5 /rm (Knutson et al. 
2012 ), 8.0 /rm (Agol et al. 2010 ), Spitzer/IRS spec¬ 
troscopy (Grillmair et al. 2008), Spitzer/IRS 16 & 
24 /im photometry (Charbonneau et al. 2008). As 
in previous studies, we consider the active trace 
gasses H 2 O, CH 4 , CO and CO 2 and otherwise as¬ 
sume a hydrogen dominated atmosphere. 

Figure 14 shows the retrieved emission models 
with retrieved abundances and log(Evidence) re¬ 
ported in table 3. Figures 15 and 16 report the 
retrieved Stage 1 & Stage 2 TP-profiles and asso¬ 
ciated Stage 1 TP-covariance, respectively. Fig¬ 
ures 17 & 18 are the marginalised and conditional 
posteriors for the active trace gasses considered. 
Figure 19 illustrates the difference between Stage 1 
and Stage 2 emission spectra. 

As described above, the Stage 1 retrieval was 
first performed using the parmetric TP-profile by 
Guillot ( 2010 ) with the additional optical opacity 
terms proposed by Line et al. ( 2012 ). Here, all 
parameters of the TP model were allowed to vary. 
We modelled the atmosphere over 25 scale heights 
sampled onto a 150 pressure-layer grid. Following 
this initial retrieval the TP-covariance was com¬ 
puted resulting in a predominantly isothermal at¬ 
mosphere down to the lowest ^20 layers (^ 0.1 
bar pressure). The obtained TP-profile is well con¬ 
strained but shows systematics (e.g. at 5 x 10 “^ 
bar) typical to this type of parametrisation. The 
posterior distributions of the active trace gases, 
figure 17 shows good constrains on abundances. 
From figure 17, we can see the retrieved values for 
methane to only constitute an upper limit. 

The Stage 2 hybrid model contained 22 free pa¬ 
rameters on its non-linearly sampled TP-profile 
grid. All parameters were allowed to vary freely 
between fully constrained {a = 1.0 in equation 16) 
and unconstrained scenarios. Figure 15 (right) 
shows the Stage 2 TP-profile with a = 0.54, indi¬ 
cating a significant shift away from the parametric 



Fig. 14.— Emission spectrum of a HD189733b, 
black circles, and Stage 1 & 2 fitting on the top 
and bottom respectively; blue: the fitted emission 
spectrum at R=1000; red squares: spectrum fit 
binned to data resolution; grey: photometric pass- 
bands for Spitzer/IRAC and MIPS. Both fits con¬ 
verged but the Stage 2 fit reached a higher global 
Evidence. 



Temperature Temperature 


Fig. 15.— TP-profiles for Stage 1 & 2 results for 
HD189733b in figure 14. Both solutions converge 
within the calculated error bar. Stage 2 features 
a significant improvement in maximum likelihood 
achieved. 
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Layers 


Fig. 16.— Correlation matrix, C, derived from 
Stage 1 TP-profile shown in figure 15. Most of 
the atmosphere is found to be best fit with an 
isothermal profile but the lowest ^ 20 atmospheric 
layers (^0.1 bar). Otherwise identical to figure 5. 



log(C'O) 


Fig. 17.— Marginalised and conditional posterior 
distributions of the trace gasses of HD189733b for 
Stage 1 fitting. All trace gases are well constrained 
but only an upper limit to methane can be found. 



\og{CO) 


Fig. 18.— Marginalised and conditional poste¬ 
rior distributions of the trace gasses of HD189733b 
for Stage 2 fitting. All trace gases are well con¬ 
strained. 


16r 


12 
' 10 


10 15 20 25 

Wavelength (/^m) 


Fig. 19.— Model fit difference between Stage 1 
and Stage 2 at spectral resolution of 1000. Differ¬ 
ences between either retrievals are of the order of 
5-8x10“^ in most wavelengths. 
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solution of Stage 1. The Stage 2 model features a 
lower tropopause pressure along with reduced wa¬ 
ter and carbon-dioxide abundances to achieve a 
very significant increase in Bayesian Evidence, see 
table 3 and discussion in section 5. The Stage 2 
TP-profile error bounds are slightly larger. Given 
a significantly higher Evidence for Stage 2^ it is 
indicative of Stage 1 TP-profile errors to be un¬ 
derestimated by the parameterisation used. The 
posteriors distributions of Stage 2 now show a con¬ 
vergence of methane beyond the upper-bound con¬ 
strained of Stage 1. 

Compared to table 3 in Line et al. (2014), we 
obtain lower abundances of CH 4 and CO 2 but 
higher values for CO. Whereas these values dif¬ 
fer from previous analyses, we find them signif¬ 
icantly closer to the chemical model predictions 
of HD189733b (e.g. Line et al. 2010; Venot et al. 
2012 ; Moses et al. 2013). Note that also Stage 1 
results exhibit the same trend of lower CH 4 and 
CO 2 abundances. We find two aspects by which 
our analysis differs from others in the literature: 
1 ) A finer and complete sampling of correlated 
likelihoods, in particular compared to maximum 
likelihood methods which are often less efficient 
with sparsely sampled data; 2) The completeness 
of molecular line-lists used. This is particularly 
true for the C/0 ratio determination using CH 4 
and CO/CO 2 as tracers. In this work we use the 
new YTlOtolO CH 4 line-list (Yurchenko & Ten¬ 
nyson 2014). 

5. Discussion 

In section 4 we have demonstrated the effi¬ 
ciency of the T-REx retrieval suite for emission 


Table 3: Retrieved abundances for hot Jupiter 
HD189733b. Top row: log Bayesian Evidence for 
Stage 1 & Stage 2 models. Alog(E) = +121.75 
indicating a very significant improvement in the 
Stage 2 fit compared to Stage 1 . 



Stage 1 Retrieval 

Stage 2 Retrieval 

log(^) 

-43.23 

78.52 

log(H20) 

-3.918 ± 0.207 

-4.978 ± 0.602 

log(CH4) 

-6.732 ± 0.719 

-6.768 ± 0.487 

log(C02) 

-3.722 ± 0.482 

-4.204 ± 0.488 

log(CO) 

-2.671 ± 1.387 

-2.689 ± 0.769 


spectroscopy using a simulated hot-Jupiter, a hot- 
SuperEarth, as well as secondary eclipse observa¬ 
tions of HD189733b, as examples. 

In all three cases the Bayesian Evidence of the 
Stage 2 retrieval is significantly higher, log(E) = 
36.20 compared to - 43.92, log(E) = 168.9 to 
75.40, log(E;) = -43.23 to 78.52 for WASP-76b, 
55 Cnc e and HD 189733b, respectively. This is 
clearly indicative of Stage 2 results being more 
robust statistically. Eollowing the adaptation of 
the Jeffrey’s scale of model evidence (Jeffreys & 
Kendall 1948) by Kass & Raftery ( 2012 ), we can 
define a strong preference for the Stage 2 model 
as Alog(E) > 5 and equally, a strong preference 
of the Stage 1 model to be Alog(E) < —5. Evi¬ 
dence differences ranging from -5 to 5 indicate a 
lesser significance. In the case of WASP-76b we 
find Alog(£’) = +80.12, Alog(E;) = +132.70 for 
55 Cnc e and Alog(E;) = +121.75 for HD189733b. 
Eurthermore the improved Stage 2 fit is illustrated 
by the better retrieval of the trace gas abundances 
(tables 1, 2 & 3). This illustrates the importance 
of accurate TP-profile retrievals and the advan¬ 
tage of a two-staged approach, especially in cases 
of low resolution and/or low S/N data. 

6. Conclusion 

In this paper we introduced the emission spec¬ 
troscopy retrieval approach for the T-REx re¬ 
trieval suite framework. Given a common code 
basis for transmission and emission retrieval, al¬ 
lows us to benefit from computational efficiencies 
and the high accuracy of molecular line-list han¬ 
dling introduced in W15. To suite the needs of the 
temperature-pressure profile retrieval, we imple¬ 
mented an extra loop unique to the emission side 
of T-REx, which allows a two-staged retrieval. We 
show that such a staged retrieval of the emission 
spectrum (and TP-profile) allows us to dynami¬ 
cally scale the complexity of the retrieval prob¬ 
lem (from a fully parameterised to a fully uncon¬ 
strained model) and has significant advantages in 
accuracy and robustness over previously available 
methods. Euture publications will extend the sen¬ 
sitivity analyses presented here to include present 
and future ground and space instrumentation and 
a wider range of planets observable in emission 
spectroscopy. 
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